home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / scheme / scm / jackal1a.lha / jacal / poly.scm < prev    next >
Encoding:
Text File  |  1992-12-23  |  15.1 KB  |  500 lines
  1. ;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
  2. ;;; Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.
  3. ;;; See the file "COPYING" for terms applying to this program.
  4.  
  5. ;;;Functions which operate on polynomials as polynomials
  6. ;;;have prefix POLY_
  7. ;;;Functions which operate on polynomials with the same major variable
  8. ;;;have prefix UNIV_
  9. ;;;Functions which operate on coefficients (without major varible)
  10. ;;;have prefix COES_
  11.  
  12. ;(proclaim '(optimize (speed 3) (compilation-speed 0)))
  13.  
  14. (define (one? n) (eqv? 1 n))
  15. (define (poly_0? p) (eqv? 0 p))
  16. (define (divides? a b) (zero? (remainder b a)))
  17.  
  18. ;;; poly is the internal workhorse data type in the form
  19. ;;; of numeric or list (var coeff0 coeff1 ...)
  20. ;;; where var is a variable and coeffn is the coefficient of var^n.
  21. ;;; coeffn is poly.  The variables currently are arranged
  22. ;;; with the reverse alphabetically with z higher order than A.
  23.  
  24. (define (poly_find-var? poly var)
  25.   (poly_find-var-if? poly (lambda (x) (eqv? var x))))
  26.  
  27. (define (poly_find-var-if? poly proc)
  28.   (cond ((number? poly) #f)
  29.     ((proc (car poly)))
  30.     (else (some (lambda (x) (poly_find-var-if? x proc)) (cdr poly)))))
  31.  
  32. ;;;POLY_VARS returns a list of all vars used in POLY
  33. (define (poly_vars poly)
  34.   (let ((elts '()))
  35.     (poly_for-each-var (lambda (v) (set! elts (adjoin v elts))) poly)
  36.     elts))
  37.  
  38. ;;;; the following functions are for internal use on the poly data type
  39.  
  40. ;;; this normalizes short polys.
  41. (define (univ_norm0 var col)
  42.   (cond ((null? col) 0)
  43.     ((null? (cdr col)) (car col))
  44.     (else (cons var col))))
  45.  
  46. (define (map-no-end-0s proc l)
  47.   (if (null? l)
  48.       l
  49.     (let ((res (proc (car l)))
  50.       (rest (map-no-end-0s proc (cdr l))))
  51.       (if (and (null? rest) (eqv? 0 res))
  52.       rest
  53.     (cons res rest)))))
  54. (define (map2c-no-end-0s proc l1 l2)
  55.   (cond ((null? l1) l2)
  56.     ((null? l2) l1)
  57.     (else
  58.      (let ((res (proc (car l1) (car l2)))
  59.       (rest (map2c-no-end-0s proc (cdr l1) (cdr l2))))
  60.        (if (and (null? rest) (eqv? 0 res))
  61.            rest
  62.          (cons res rest))))))
  63. (define (univ-w/o-lt p)
  64.   (univ_norm0 (car p) (map-no-end-0s identity (butlast (cdr p) 1))))
  65. (define (make-list-length lst len fill)
  66.   (let ((ld (- len (length lst))))
  67.     (cond ((negative? ld) (butlast lst (- ld)))
  68.       (else (append lst (make-list ld fill))))))
  69. ;;; replaces each var in poly with (proc var)
  70. (define (poly_do-vars proc poly)
  71.   (if (number? poly)
  72.       poly
  73.     (cons (proc (car poly))
  74.       (map (lambda (b) (poly_do-vars proc b))
  75.            (cdr poly)))))
  76. (define (polys_do-vars proc polys)
  77.   (bunch_map (lambda (poly) (poly_do-vars proc poly)) polys))
  78. (define (poleqns_do-vars proc licits)
  79.   (bunch_map (lambda (poly) (poly_do-vars proc (licit->poleqn poly))) licits))
  80. ;;; This can call proc more than once per var
  81. (define (poly_for-each-var proc poly)
  82.   (cond ((number? poly))
  83.     (else
  84.      (proc (car poly))
  85.      (for-each (lambda (b) (poly_for-each-var proc b))
  86.            (cdr poly)))))
  87. (define (polys_for-each-var proc polys)
  88.   (bunch_for-each (lambda (poly) (poly_for-each-var proc poly)) polys))
  89. (define (ipow-by-squaring x n acc proc)
  90.   (cond ((zero? n) acc)
  91.     ((one? n) (proc acc x))
  92.     (else (ipow-by-squaring (proc x x)
  93.                 (quotient n 2)
  94.                 (if (even? n) acc (proc acc x))
  95.                 proc))))
  96.  
  97. (define (poly_add-const term p2)
  98.   (cons (car p2) (cons (poly_+ term (cadr p2)) (cddr p2))))
  99. (define (poly_+ p1 p2)
  100.   (cond ((and (number? p1) (number? p2)) (+ p1 p2))
  101.     ((number? p1) (if (zero? p1) p2 (poly_add-const p1 p2)))
  102.     ((number? p2) (if (zero? p2) p1 (poly_add-const p2 p1)))
  103.     ((eq? (car p1) (car p2))
  104.      (univ_norm0 (car p1) (map2c-no-end-0s poly_+ (cdr p1) (cdr p2))))
  105.     ((var_> (car p2) (car p1)) (poly_add-const p1 p2))
  106.     (else (poly_add-const p2 p1))))
  107.  
  108. (define (univ_* p1 p2)
  109.   (let ((res (make-list (+ (length (cdr p1)) (length (cdr p2)) -1) 0)))
  110.     (do ((rpl res (cdr rpl))
  111.      (a (cdr p1) (cdr a)))
  112.     ((null? a) (cons (car p1) res))
  113.     (do ((b (cdr p2) (cdr b))
  114.          (rp rpl (cdr rp)))
  115.         ((null? b))
  116.         (set-car! rp (poly_+ (poly_* (car a) (car b)) (car rp)))))))
  117. (define (poly_times-const term p2)
  118.   (cons (car p2) (map (lambda (x) (poly_* term x)) (cdr p2))))
  119. (define (poly_* p1 p2)
  120.   (cond ((and (number? p1) (number? p2)) (* p1 p2))
  121.     ((number? p1) (cond ((zero? p1) 0)
  122.                 ((one? p1) p2)
  123.                 (else (poly_times-const p1 p2))))
  124.     ((number? p2) (cond ((zero? p2) 0)
  125.                 ((one? p2) p1)
  126.                 (else (poly_times-const p2 p1))))
  127.     ((eq? (car p1) (car p2)) (univ_* p1 p2))
  128.     ((var_> (car p2) (car p1)) (poly_times-const p1 p2))
  129.     (else (poly_times-const p2 p1))))
  130.  
  131. (define (poly_negate p) (poly_* -1 p))
  132. (define (poly_- p1 p2) (poly_+ p1 (poly_negate p2)))
  133.  
  134. (define (univ_/? u v)
  135.   (let ((r (list->vector (cdr u)))
  136.     (m (length (cddr u)))
  137.     (n (length (cddr v)))
  138.     (vn (car (last-pair v)))
  139.     (q '()))
  140.     (do ((k (- m n) (- k 1))
  141.      (qk (poly_/? (vector-ref r m) vn)
  142.          (and (> k 0) (poly_/? (vector-ref r (+ n k -1)) vn))))
  143.     ((not qk)
  144.      (and (< k 0)
  145.           (do ((k (- n 2) (- k 1)))
  146.           ((or (< k 0) (not (poly_0? (vector-ref r k))))
  147.            (< k 0)))
  148.           (univ_norm0 (car u) q)))
  149.       (set! q (cons qk q))
  150.       (let ((qk- (poly_negate qk)))
  151.     (do ((j (+ n k -1) (- j 1)))
  152.         ((< j k))
  153.       (vector-set! r j (poly_+ (vector-ref r j)
  154.                    (poly_* (list-ref v (+ (- j k) 1)) qk-))))))))
  155.  
  156. ;;; POLY_/? returns U / V if V divides U, otherwise returns #f
  157. (define (poly_/? u v)
  158.   (cond ((equal? u v) 1)
  159.     ((eqv? 0 u) 0)
  160.     ((eqv? 0 v) #f)
  161.     ((unit? v) (poly_* v u))
  162.     ((and (number? u) (number? v)) (and (divides? v u) (quotient u v)))
  163.     ((number? v) (univ_/? u (const_promote (car u) v)))
  164.     ((number? u) #f)
  165.     ((eq? (car u) (car v)) (univ_/? u v))
  166.     ((var_> (car u) (car v))
  167.      (univ_/? u (const_promote (car u) v)))
  168.     (else #f)))
  169.  
  170. (define (univ_/ dividend divisor)
  171.   (or (univ_/? dividend divisor)
  172.       (math-error divisor " does not udivide " dividend)))
  173.  
  174. (define (poly_/ dividend divisor)
  175.   (or (poly_/? dividend divisor)
  176.       (math-error divisor " does not divide " dividend)))
  177.  
  178. (define (univ_monomial coeff n var)
  179.   (cond ((eq? 0 coeff) 0)
  180.     ((>= 0 n) coeff)
  181.     (else
  182.      (cons var
  183.            ((lambda (x) (set-car! (last-pair x) coeff) x)
  184.         (make-list (+ 1 n) 0))))))
  185.  
  186. (define (poly_degree p var)
  187.    (cond ((number? p) 0)
  188.      ((eq? var (car p)) (length (cddr p)))
  189.      ((var_> var (car p)) 0)
  190.      (else (reduce-init (lambda (m c) (max m (poly_degree c var)))
  191.                 0
  192.                 (cdr p)))))
  193.  
  194. ;(define (leading-coeff p var) (poly_coeff p var (poly_degree p var)))
  195. ;(define (monic? u var) (one? (leading-coeff u var)))
  196.  
  197. (define (poly_^ x n)
  198.   (if (number? x)
  199. ;      (expt x n)
  200.     (ipow-by-squaring x n 1 *)
  201.     (ipow-by-squaring x n 1 poly_*)))
  202.  
  203. ;;;; Routines used in normalizing IMPL polynomials
  204.  
  205. (define (leading-number p)
  206.   (if (number? p) p (leading-number (car (last-pair p)))))
  207.  
  208. ;;; This canonicalizes polys with respect to sign by forcing the
  209. ;;; numerical coefficient of the a certain term to always be positive.
  210. (define (signcan p)
  211.   (if (> (leading-number p) 0) (poly_negate p) p))
  212.  
  213. (define (shorter? x y) (< (length x) (length y)))
  214.  
  215. (define (univ_degree p var)
  216.   (if (or (number? p) (not (eq? (car p) var))) 0 (length (cddr p))))
  217.  
  218. ;;; THE NEXT 3 ROUTINES FOR SUBRESULTANT GCD ASSUME THAT THE ARGUMENTS
  219. ;;; ARE POLYNOMIALS WITH THE SAME MAJOR VARIABLE.
  220. ;;;  THESE TWO ROUTINES ASSUME THAT THE FIRST ARGUMENT IS OF GREATER
  221. ;;;  OR EQUAL ORDER THAN THE SECOND.
  222. ;;; These algorithms taken from:
  223. ;;; Knuth, D. E.,
  224. ;;; The Art Of Computer Programming, Vol. 2: Seminumerical Algorithms,
  225. ;;; Addison Wesley, Reading, MA 1969.
  226. ;;; Pseudo Remainder
  227.  
  228. ;;; This returns a list of the pseudo quotient and pseudo remainder.
  229. (define (univ_pdiv u v)
  230.   (let* ((r (list->vector (cdr u)))
  231.      (m (length (cddr u)))
  232.      (n (length (cddr v)))
  233.      (vn (car (last-pair v)))
  234.      (q (make-vector (+ (- m n) 1) 1)))
  235.     (do ((tt (- (- m n) 1) (- tt 1))
  236.      (k 1 (+ 1 k))
  237.      (vnp 1))
  238.     ((< tt 0))
  239.     (set! vnp (poly_* vnp vn))
  240.     (vector-set! q k vnp)
  241.     (vector-set! r tt (poly_* (vector-ref r tt) vnp)))
  242.     (do ((k (- m n) (- k 1))
  243.      (rnk 0))
  244.     ((< k 0))
  245.       (set! rnk (poly_negate (vector-ref r (+ n k))))
  246.       (do ((j (+ n k -1) (- j 1)))
  247.       ((< j k))
  248.     (vector-set! r j (poly_+ (poly_* (vector-ref r j) vn)
  249.                  (poly_* (list-ref v (+ (- j k) 1)) rnk)))))
  250.     (list
  251.      (do ((k (- m n) (+ -1 k))
  252.       (end '() (cons (poly_* (vector-ref r (+ n k))
  253.                  (vector-ref q k)) end)))
  254.      ((zero? k) (univ_norm0 (car u) (cons (vector-ref r n) end))))
  255.      (do ((j (- n 1) (- j 1))
  256.       (end '()))
  257.      ((< j 0) (univ_norm0 (car u) end))
  258.      (if (not (and (null? end) (eqv? 0 (vector-ref r j))))
  259.          (set! end (cons (vector-ref r j) end)))))))
  260. (define (poly_pdiv dividend divisor var)
  261.   (let ((pd1 (poly_degree dividend var))
  262.     (pd2 (poly_degree divisor var)))
  263.     (cond ((< pd1 pd2) (list 0 dividend))
  264.       ((zero? (+ pd1 pd2))
  265.        (list (quotient dividend divisor) (remainder dividend divisor)))
  266.       ((zero? pd1) (list 0 dividend))
  267.       ((zero? pd2)
  268. ;;; This should work but doesn't.
  269. ;;;       (map univ_demote (univ_pdiv (promote var dividend)
  270. ;;;                       (const_promote var divisor)))
  271.        (list 0 dividend))
  272.       (else
  273.        (map univ_demote (univ_pdiv (promote var dividend)
  274.                        (promote var divisor)))))))
  275. (define (univ_prem u v)
  276.   (let* ((r (list->vector (cdr u)))
  277.      (m (length (cddr u)))
  278.      (n (length (cddr v)))
  279.      (vn (car (last-pair v))))
  280.     (do ((k (- (- m n) 1) (- k 1))
  281.      (vnp 1))
  282.     ((< k 0))
  283.       (set! vnp (poly_* vnp vn))
  284.       (vector-set! r k (poly_* (vector-ref r k) vnp)))
  285.     (do ((k (- m n) (- k 1))
  286.      (rnk 0))
  287.     ((< k 0))
  288.       (set! rnk (poly_negate (vector-ref r (+ n k))))
  289.       (do ((j (+ n k -1) (- j 1)))
  290.       ((< j k))
  291.     (vector-set! r j (poly_+ (poly_* (vector-ref r j) vn)
  292.                  (poly_* (list-ref v (+ (- j k) 1)) rnk)))))
  293.     (do ((j (- n 1) (- j 1))
  294.      (end '()))
  295.     ((< j 0) (univ_norm0 (car u) end))
  296.       (if (and (null? end) (eqv? 0 (vector-ref r j)))
  297.       #f
  298.       (set! end (cons (vector-ref r j) end))))))
  299.  
  300. ;;; Pseudo Remainder Sequence
  301. (define (univ_prs u v)
  302.   (let ((var (car u))
  303.     (g 1)
  304.     (h 1)
  305.     (delta 0))
  306.     (do ((r (univ_prem u v) (univ_prem u v)))
  307.     ((eqv? 0 (univ_degree r var))
  308.      (if (eqv? 0 r) v r))
  309.       (set! delta (- (univ_degree u var) (univ_degree v var)))
  310.       (set! u v)
  311.       (set! v (univ_/ r (const_promote (car r) (poly_* g (poly_^ h delta)))))
  312.       (set! g (car (last-pair u)))
  313.       (set! h (cond ((one? delta) g)
  314.             ((zero? delta) h)
  315.             (else (poly_/ (poly_^ g delta)
  316.                   (poly_^ h (+ -1 delta)))))))))
  317.  
  318. (define (univ_gcd u v)
  319.   (let* ((cu (univ_cont u))
  320.      (cv (univ_cont v))
  321.      (c (poly_gcd cu cv))
  322.      (ppu (poly_/ u cu))
  323.      (ppv (poly_/ v cv))
  324.      (ans (if (shorter? ppv ppu)
  325.           (univ_prs ppu ppv)
  326.         (univ_prs ppv ppu))))
  327.     (if (zero? (univ_degree ans (car u)))
  328.     c
  329.       (poly_* c (univ_primpart ans)))))
  330.  
  331. (define (poly_gcd p1 p2)
  332.   (cond ((equal? p1 p2) p1)
  333.     ((and (number? p1) (number? p2)) (gcd p1 p2))
  334.     ((number? p1) (if (zero? p1)
  335.               p2
  336.             (apply poly_gcd* p1 (cdr p2))))
  337.     ((number? p2) (if (zero? p2)
  338.               p1
  339.             (apply poly_gcd* p2 (cdr p1))))
  340.     ((eq? (car p1) (car p2)) (univ_gcd p1 p2))
  341.     ((var_> (car p2) (car p1)) (apply poly_gcd* p1 (cdr p2)))
  342.     (else (apply poly_gcd* p2 (cdr p1)))))
  343.  
  344. (define (poly_gcd* . li)
  345.   (let ((nums (remove-if-not number? li)))
  346.     (if (null? nums)
  347.     (reduce poly_gcd li)
  348.     (let ((gnum (reduce gcd nums)))
  349.       (if (= 1 gnum) 1
  350.           (reduce-init poly_gcd gnum (remove-if number? li)))))))
  351.  
  352. (define (univ_cont p) (apply poly_gcd* (cdr p)))
  353. (define (univ_primpart p) (poly_/ p (univ_cont p)))
  354. (define (poly_num-cont p)
  355.   (if (number? p) p
  356.       (do ((l (cdr p) (cdr l))
  357.        (n (poly_num-cont (cadr p))
  358.           (gcd n (poly_num-cont (cadr l)))))
  359.       ((or (= 1 n) (null? (cdr l))) n))))
  360.  
  361. (define (list-ref? l n)
  362.   (cond ((null? l) #f)
  363.     ((zero? n) (car l))
  364.     (else (list-ref? (cdr l) (- n 1)))))
  365.  
  366. (define (univ_coeff p ord) (or (list-ref? (cdr p) ord) 0))
  367. (define (poly_coeff p var ord)
  368.   (cond ((or (number? p) (var_> var (car p)))
  369.      (if (zero? ord) p 0))
  370.     ((eq? var (car p)) (univ_coeff p ord))
  371.     (else
  372.      (univ_norm0 (car p)
  373.              (map-no-end-0s (lambda (c) (poly_coeff c var ord))
  374.                     (cdr p))))))
  375.  
  376. (define (poly_subst0 old e) (poly_coeff e old 0))
  377.  
  378. (define const_promote list)
  379.  
  380. (define (promote var p)
  381.   (if (eq? var (car p))
  382.       p
  383.       (let ((dgr (poly_degree p var)))
  384.     (do ((i dgr (+ -1 i))
  385.          (ol (list (poly_coeff p var dgr))
  386.          (cons (poly_coeff p var (+ -1 i)) ol)))
  387.         ((zero? i) (cons var ol))))))
  388.  
  389. ;;;this is bummed if v has higher priority than any variable in (cdr p)
  390. (define (univ_demote p)
  391.   (if (number? p)
  392.       p
  393.     (let ((v (car p)))
  394.       (if (every (lambda (cof) (or (number? cof) (var_> v (car cof))))
  395.          (cdr p))
  396.       p
  397.     (poly_+ (cadr p)
  398.         (do ((trms (cddr p) (cdr trms))
  399.              (sum 0)
  400.              (mon (list v 0 1) (cons v (cons 0 (cdr mon)))))
  401.             ((null? trms) sum)
  402.             (set! sum (poly_+ sum (poly_* mon (car trms))))))))))
  403.  
  404. (define (sylvester p1 p2 var)
  405.   (set! p1 (promote var p1))
  406.   (set! p2 (promote var p2))
  407.   (let ((d1 (univ_degree p1 var))
  408.     (d2 (univ_degree p2 var))
  409.     (m (list)))
  410.     (do ((i d1 (+ -1 i))
  411.      (row (nconc (make-list (+ -1 d1) 0) (reverse (cdr p2)))
  412.           (append (cdr row) (list 0))))
  413.     ((<= i 1) (set! m (cons row m)))
  414.     (set! m (cons row m)))
  415.     (do ((i d2 (+ -1 i))
  416.      (row (nconc (make-list (+ -1 d2) 0) (reverse (cdr p1)))
  417.           (append (cdr row) (list 0))))
  418.     ((<= i 1) (set! m (cons row m)))
  419.     (set! m (cons row m)))
  420.     m))
  421.  
  422. ;;; Bareiss's integer preserving gaussian elimination.
  423. ;;; Bareiss, E.H.: Sylvester's identity and multistep
  424. ;;; integer-preserving Gaussian elimination. Mathematics of
  425. ;;; Computation 22, 565-578, 1968.
  426. ;;; as related by:
  427. ;;; Akritas, A.G.: Exact Algorithms for the Matrix-Triangulation
  428. ;;; Subresultant PRS Method.  Computers and Mathematics, 145-155.
  429. ;;; Springer Verlag, 1989.
  430. (define (bareiss m)
  431.   4)
  432.  
  433. (define (poly_resultant p1 p2 var)
  434.   (let ((u1 (promote var p1))
  435.     (u2 (promote var p2)))
  436.     (or (not (zero? (univ_degree u1 var)))
  437.     (not (zero? (univ_degree u2 var)))
  438.     (math-error var " does not appear in " p1 " or " p2))
  439.     (let ((res (cond ((zero? (univ_degree u1 var)) p1)
  440.              ((zero? (univ_degree u2 var)) p2)
  441.              ((shorter? u1 u2) (univ_prs u2 u1))
  442.              (else (univ_prs u1 u2)))))
  443.       (if (zero? (univ_degree res var)) res
  444.       0))))
  445.  
  446. (define (poly_mod_number poly modulus)
  447.   (if (number? poly)
  448.       (modulo poly modulus)
  449.     (cons (car poly)
  450.       (map-no-end-0s (lambda (x) (poly_mod_number x modulus))
  451.              (cdr poly)))))
  452.  
  453. (define (poly_prem dividend divisor)
  454.   (if (number? divisor)
  455.       (poly_mod_number dividend divisor)
  456.       (let ((var (car divisor)))
  457.     (if (> (poly_degree divisor var) (poly_degree dividend var))
  458.         dividend
  459.         (univ_demote (univ_prem (promote var dividend)
  460.                     (promote var divisor)))))))
  461.  
  462. ;;;; VERIFICATION TESTS
  463. (define (poly_test)
  464.   (define a (sexp->var 'A))
  465.   (define b (sexp->var 'B))
  466.   (define c (sexp->var 'C))
  467.   (define x (sexp->var 'X))
  468.   (define y (sexp->var 'Y))
  469.   (test (list a 0 -2)
  470.     poly_gcd
  471.     (list a 0 -2)
  472.     (list a 0 0 -2))
  473.   (test (list x (list a 0 1) 1)
  474.     poly_gcd
  475.     (list x (list a 0 0 -1) 0 1)
  476.     (list x (list a 0 0 1) (list a 0 2) 1))
  477.   (test (list x 0 (list a 0 1))
  478.     poly_gcd
  479.     (list x 0 (list a 0 0 1))
  480.     (list x 0 0 (list a 0 1)))
  481.   (test (list x (list b 0 0 1) 0 (list b 1 2) (list a 0 1) 1)
  482.     poly_resultant
  483.     (list y (list x (list b 0 1) 0 1) (list x 0 1))
  484.     (list y (list x 1 (list a 0 1)) 0 1)
  485.     y)
  486.   (test (list y (list b 0 0 1) 0 (list b 1 2) (list a 0 1) 1)
  487.     poly_resultant
  488.     (list y (list b 0 1) (list x 0 1) 1)
  489.     (list y (list x 1 0 1) (list a 0 1))
  490.     x)
  491.   '(test 1
  492.      poly_gcd
  493.      (list x -5 2 8 -3 -3 1 1)
  494.      (list x 21 -9 -4 5 3))
  495.   '(test 1
  496.      poly_gcd
  497.      (list x -5 2 8 -3 -3 0 1 0 1)
  498.      (list x 21 -9 -4 0 5 0 3))
  499.   'done)
  500.